Anxiety, concerns and COVID-19: Cross-country perspectives from families and individuals with neurodevelopmental conditions

Background The COVID-19 pandemic had a major impact on the mental health and well-being of children with neurodevelopmental conditions (NDCs) and of their families worldwide. However, there is insufficient evidence to understand how different factors (e.g., individual, family, country, children) have impacted on anxiety levels of families and their children with NDCs developed over time. Methods We used data from a global survey assessing the experience of 8043 families and their children with NDCs (mean of age (m) = 13.18 years, 37% female) and their typically developing siblings (m = 12.9 years, 45% female) in combination with data from the European Centre for Disease Prevention and Control, the University of Oxford, and the Central Intelligence Agency (CIA) World Factbook, to create a multilevel data set. Using stepwise multilevel modelling, we generated child-, family- and country-related factors that may have contributed to the anxiety levels of children with NDCs, their siblings if they had any, and their parents. All data were reported by parents. Results Our results suggest that parental anxiety was best explained by family-related factors such as concerns about COVID-19 and illness. Children’s anxiety was best explained by child-related factors such as children’s concerns about loss of routine, family conflict, and safety in general, as well as concerns about COVID-19. In addition, anxiety levels were linked to the presence of pre-existing anxiety conditions for both children with NDCs and their parents. Conclusions The present study shows that across the globe there was a raise in anxiety levels for both parents and their children with NDCs because of COVID-19 and that country-level factors had little or no impact on explaining differences in this increase, once family and child factors were considered. Our findings also highlight that certain groups of children with NDCs were at higher risk for anxiety than others and had specific concerns. Together, these results show that anxiety of families and their children with NDCs during the COVID-19 pandemic were predicted by very specific concerns and worries which inform the development of future toolkits and policy. Future studies should investigate how country factors can play a protective role during future crises.

. Sample sizes by country. Number of unique families, children, and diagnosis groups.

B.1 Data analysis plan
Analysis of the data consisted of two general steps, (1) pre-processing, and (2) stepwise multilevel regressions of parent and child anxiety. These two general steps consisted of further specific steps, as summarized visually in Figure B1. Figure B1. Pipeline of data analysis, consisting of (1) pre-processing followed by (2) stepwise multilevel regressions for parent and child anxiety. IV = independent variables.

B.2 Pre-processing
First, all country samples were merged and cleaned for language, notation, and formatting consistency. Second, respondents and variables that were unusable for analysis were removed a priori. This removal included (a) respondents that disagreed to have their data analyzed, (b) respondents that completed less than 25% of the questionnaire, (c) respondents that failed a check that separated human users from automated bots, (d) respondents for which the primary diagnosis of their NDCs child was missing, (e) variables that contained open text responses or raw dates, and (f) variables that had over 70% of their observations missing. This exclusion reduced the number of families from 22,741 to 8,043, the number of countries from 128 to 78, and the number of variables from 406 to 268. Third, country variables were added to the cleaned data that measured national COVID19 evolution, national government responses to the pandemic, and structural country descriptors (e.g., size, population). COVID19 evolution statistics during the time of the study (until 31.07.2020) were sourced from the European Center for Disease Prevention and Control. 1 From these statistics, we calculated for each country total deaths, total deaths per million inhabitants, peak deaths, peak deaths per million, date of the first death, date of peak deaths per million. Peak deaths and peak deaths per million were calculated from a 7-day moving average of deaths per million, to reduce the impact of outlying values. From the calculated evolution statistics, two new variables were derived to express for each family their distance to the international pandemic start, and their distance to their national pandemic peak. Distance to the international pandemic start was calculated as the difference in days between 11.03.2020 2 and the day the family completed the questionnaire (i.e., the "now" moment). Distance to the national pandemic peak was calculated as the difference in days between the day peak deaths were achieved nationally and the day the family completed the questionnaire. 3 Multilevel regressions subsequently controlled for these two distance variables, to adjust for the fact that different families completed the questionnaire at different times during the pandemic.
Government response statistics were sourced from the Oxford COVID-19 Government Response Tracker (Hale et al., 2021 4 ), which included the tracked government response to the pandemic in multiple areas (e.g., public transport, teleworking, financial aid) and four aggregated indexes as general measures of response severity (including the stringency index -containment and closure policies, sometimes referred to as lockdown policies). Structural country descriptors were sourced from the CIA World Factbook, 5 and included variables on geographic, demographic, economic and political aspects of each country.
Fourth, a univariate outlier screening was conducted on the remaining data to check for impossible or unusual values. For categorical variables, levels with low occurrence were identified as outlying. For continuous variables, boxplots were used to identify observations beyond 1.5 times the interquartile distance of the rank-ordered data. Impossible values were converted to missing, while outlying values were simply flagged for later analyses.
Fifth, the pre-processed and expanded data were restructured from "wide format" to "long format" for compatibility with multilevel regression. In long format, each respondent/family had between 3 or 6 rows of repeated measures data, corresponding to the child (NDCs and TD) × time (before, start, now) repeated measures structure. As such, the full hierarchical structure of the data consisted of time points (before, start, now) nested in children (NDCs and TD), children nested in families, and families nested in countries. The presence of this data hierarchy motivated the use of multilevel models for stepwise analysis.

B.3 Stepwise multilevel regression
The final long-format data were entered into two stepwise multilevel regressions, with as dependent variables (DV) parent anxiety and child anxiety, respectively. For both DVs, modelling proceeded identically as a blocked stepwise regression search. Four blocks of independent variables (IVs) were added incrementally to obtain a final model, (1) base variables, (2) country variables, (3) family variables, and (4) child variables. The base model included the primary design IVs, which were the diagnosis group variable (ASD = Autism Spectrum Disorder, ADHD = Attention-Deficit/Hyperactivity Disorder, DLD = Developmental Language Disorder, DS = Down syndrome, WS = Williams syndrome, ID = Intellectual Disability-Not otherwise specified, TDS = Typically developing sibling) interacting with categorical time ("before", "start", "now"). As well, the base model included controlling covariates for demographical parent and child aspects (e.g., age, gender), each family's distance to the international pandemic start, and the distance to their national pandemic peak. Subsequent IV blocks add further predictors of anxiety going from the macro-level of influence (i.e., country aspects) to the microlevel of influence (i.e., child aspects). The full list of variables for blocks 1-4 can be found in Table 2 in the main paper.
For each block, fitting the multilevel model consisted of four steps, (a) fixed effects selection, (b) random effects selection, (c) model diagnostics and inspection, and (d) ANOVA breakdown of effects. However, step (a) was omitted for the base model, due to these IVs representing structural information that we wished to control for, no matter its statistical significance.
Fixed effects selection. Within blocks 2-4, important fixed effects were identified by a forward selection pass, followed by backward elimination pass: • Forward selection pass: At each step of the pass, the IV entered the model with the lowest p-value that surpassed 0.005 in significance. This was repeated until no new IV reached the threshold for inclusion. • Backward elimination pass: At each step of the pass, the IV in the model with the highest p-value that was less significant than 0.05 was removed. This was repeated until no further IV in the model reached the threshold for elimination.
For the final set of selected IVs within each block, we further investigated possible two-way interactions between the selected IVs and either time or diagnosis group.
Random effects selection. Once the important fixed effects had been identified in a block, important random effects were identified in that block. Random effects in a multilevel model reflect different sources of sampling error (e.g., repeated sampling of populations such as countries and families) and enabled us to consider the hierarchical nesting in the data. For block 1, we considered a base random effects structure consisting of a random intercept for countries and a random intercept for families, as well as random slopes for categorical time within those two intercepts. This base structure accounted for the possibility that there were baseline group and time differences due to random sampling of families and countries and was not subjected to further reduction. For blocks 2-4, important random slopes were identified with a stepwise forward and backward search, identical to the strategy conducted for the fixed effects. This time, Akaike's Information Criterion (AIC) was used as the criterion to add or remove random effects, with a relative increase or decrease in more than 2 AIC points considered as the threshold. AIC selection was preferred for the random effects due to complications related to inferential tests on variance parameters (see Fitzmaurice, Laird & Ware, 2004).
Model diagnostics. The model that was retained for each block by the fixed and random effects selection was first checked for its stability and statistical assumptions by inspecting (a) collinearity diagnostics, (b) influence diagnostics, and (c) residual diagnostics (linearity, normality, and homoscedasticity). For collinearity, we inspected variance inflation factors (VIF) for fixed effects, removing any IVs with a VIF that exceeded 10 (Kutner et al., 2005). For influence diagnostics, we inspected boxplots of model residuals and random effects, the former for identifying outlying cases, the latter for identifying outlying levels of the data hierarchy (e.g., outlying countries). For residual diagnostics, we first computed decorrelated residuals according to the procedure recommended by Fitzmaurice, Laird, and Ware (2004). Then, we plotted these transformed residuals against transformed fitted values to check for issues with non-normality, non-constant variance, and non-linear associations.
ANOVA breakdown. Once a model had been diagnosed to be stable in a block, a Type II ANOVA breakdown using F-tests was calculated to inspect statistically significant effects. For multi-parameter F-tests (e.g., diagnosis group effect), pairwise contrasts were calculated using t-tests. All inferential tests were conducted at a reduced significance level of α = 0.001, following recent recommendations to reduce the rate of false positives in social sciences research (Benjamin et al., 2018). Degrees of freedom for inferential F-and t-tests were corrected for the presence of random effects according to the method of Satterthwaite (Fitzmaurice et al., 2004), which yields appropriate fractional degrees of freedom.
As a measure for goodness-of-fit, we calculated marginal R 2 for each added block of IVs. This allowed us to quantify cumulative proportion of (fixed) variance explained by each incrementally added block of IVs (i.e., country, family, child). As a measure of effect size, we calculated standardized regression coefficients for all effects. 6 For effects with more than two categorical levels, we selected the coefficient that captured the largest standardized difference observed between two levels (Maxwell & Delaney, 2004). Finally, effect plots were generated to visualize relevant interactions and main effects.